/*
 Copyright 2014-Present Algorithms in Motion LLC
 
 This file is part of FDTD++.
 
 FDTD++ is proprietary software: you can use it and/or modify it
 under the terms of the Algorithms in Motion License as published by
 Algorithms in Motion LLC, either version 1 of the License, or (at your
 option) any later version.
 
 FDTD++ is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 Algorithms in Motion License for more details.
 
 You should have received a copy of the Algorithms in Motion License
 along with FDTD++. If not, see <http://www.aimotionllc.com/licenses/>.
*/
// CREATED    : 10/11/2015
// LAST UPDATE: 10/12/2015

#include "statsxx/machine_learning/restricted_Boltzmann_machine/RBM.hpp"

// jScience
#include "jScience/linalg.hpp" // Vector<>, transpose()

// stats++
#include "statsxx/distribution.hpp" // distribution::binomial()
#include "statsxx/machine_learning/activation_functions.hpp" // activation_function::Logistic


// basic algorithm for one step of Gibbs sampling, starting from the visible layer
//
// note: the input v is taken as provided (possibly continuous); upon return, is it sampled if mean_field == false
//
inline Vector<double> restricted_Boltzmann_machine::RBM::Gibbs_vhv(
                                                                   Vector<double> v,          // initial visible data
                                                                   const bool mean_field      // mean-field approximation
                                                                   ) const
{
    activation_function::Logistic logistic;
    
    // sample h
    // note: [nh x nv] * [nv x 1] = [nh x 1]
    Vector<double> h = logistic.f(this->c + this->W*v);
    
    // note: the hidden layer is ALWAYS sampled
    h = distribution::binomial<double>(1, h);
    
    // sample v
    // note: [nh x nv]^T * [nh x 1] = [nv x 1]
    v = logistic.f(this->b + transpose(this->W)*h);
    
    if(!mean_field)
    {
        v = distribution::binomial<double>(1, v);
    }

    return v;
}


